function r = simtnorm_double(a,b,u)

% Left and right limits
xmin = -5;
xmax =  5;

xminmin = -10;
xmaxmax = 10;

if size(u, 1)~=size(a, 1)
    error('must be a column vector');
end

% simulate RV that partially overlaps with [xmin, xmax]
r = nan(size(u));
indmin = (a>xmin & a<xmax);

if any(indmin)
    amin = a(indmin);    
    bmin = min(xmax, b(indmin));    
    r(indmin, :) = norminv(normcdf(amin) +...
        u(indmin, : ).*(normcdf(bmin) - normcdf(amin)));    
end

indtail = (a<=xmin & b>xmax);
if any(indtail)
    amax = max(xmin, a(indtail));
    bmin = min(xmax, b(indtail)); 
    r(indtail, :) = norminv(normcdf(amax) +...
        u(indtail, :).*(normcdf(bmin) - normcdf(amax)));
end

indmax = (a<=xmin & b<xmax & b>xmin);
if any(indmax)
    amax = max(xmin, a(indmax));
    bmax = b(indmax);
    r(indmax, :) = norminv(normcdf(amax) +...
        u(indmax, :).*(normcdf(bmax) - normcdf(amax)));
end

indn = ~indmin & ~indmax & ~indtail & a<0;
nn = sum(indn);
indp = ~indmin & ~indmax & ~indtail & a>0;

if any([indn; indp])
    abtail = [[abs(b(indn)), abs(b(indn))+5]; [a(indp), a(indp)+5]];
    un = 1-u(indn, :);
    utail = [un; u(indp, :)];
    q = @(x)normcdf(-x)./normpdf(x);
    qa = q(abtail(:, 1));
    qb = q(abtail(:, 2));
    c = qa.*(1-utail) + qb.*utail.*exp((abtail(:, 1).^2-abtail(:, 2).^2)/2);
    dx = 1000;
    z = 1-utail+utail.*exp((abtail(:, 1).^2-abtail(:, 2).^2)/2);
    x = sqrt(abtail(:, 1).^2-2*log(z));

    ct = 0;
    while dx>1e-10 && ct<20
        z = z-x.*(z.*q(x)-c);
        xnew = sqrt(abtail(:, 1).^2-2*log(z));        
        dx = max(abs(xnew(:)-x(:))./x(:));        
        x = xnew;        
        ct = ct+1;
    end
    
    if any(indn)
        r(indn, :) = -x(1:nn, :);
    end
    
    if any(indp)
        r(indp, :) = x(nn+1:end, :);
    end
end

indminmin = (a<=xminmin & b<=xminmin);
indmaxmax = (a>=xmaxmax & b>=xmaxmax);
if any(indminmin)
    r(indminmin, :) = b(indminmin) - repmat(abs(b(indminmin))/10, [1 size(u, 2)]);
end
if any(indmaxmax)    
    r(indmaxmax, :) = a(indmaxmax) + repmat(abs(a(indmaxmax))/10, [1 size(u, 2)]);
end























